--- title: Buffon's needles author: Ty Tuff date: '2018-12-07' slug: pie-from-matches categories: [] tags: [] ---

Make a match

match_stick_length <- 10
number_of_matches_in_a_box <- 1000
distance_between_surface_lines <- 20

point_1 <- c(-(match_stick_length/2),(match_stick_length/2))
point_2 <- c(0.01,0.01)
single_match <-cbind(point_1, point_2)
rownames(single_match)<- c("point_1","point_2")
colnames(single_match)<- c("x","y")
single_match
##          x    y
## point_1 -5 0.01
## point_2  5 0.01
library(rgl)
library(shiny)


plotids <- with(iris, 
               
      readSTL( "~/Downloads/STL/Bolt Big.STL", plot=TRUE, col="grey")

                )
rglwidget(elementId = "plot3drgl222")
plotids <- with(iris, 
               plot3d(x=c(-5,5), y=c(0,0), z=c(0,0), box=FALSE, axes=FALSE, type="n")
                )
plotids <- with(iris, 
               lines3d(x=c(-5,5), y=c(0,0), z=c(0,0), col="black", lwd=5)
                )
plotids <- with(iris, 
               spheres3d(x=c(-5,5), y=c(0,0), z=c(0,0), col="green", radius=1)
                )
rglwidget(elementId = "plot3drgl3")

Throw the match

First rotate around the origin

#install.packages("DescTools")
library(DescTools)
single_match_coords <- xy.coords(single_match)
after_rotation <-Rotate(single_match_coords$x, single_match_coords$y, 90*pi/180)


point_coords = single_match

rotate_needle <- function(angle=45, point_coords){
  angle_input <- angle*pi/180
xy <- as.matrix(point_coords)
    cos.angle <- cos(angle_input)
    sin.angle <- sin(angle_input)
    after_rotation <- point_coords
    after_rotation[1:2,1:2] <- rep(NA,4)
    after_rotation[1,] <- t(as.matrix(xy[1,])) %*% t(matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2))
    after_rotation[2,] <- t(as.matrix(xy[2,])) %*% t(matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2))
    return(after_rotation)
}

after_rotation <- rotate_needle(15, single_match)
polar_circle <- matrix(NA,4,361)

for (k in seq(0,361,length.out=360)){  ## long, angle measured counter clockwise looking down the north pole        #from the x-axis. 
    
    #i <- 89
    #k <- 55
    
    
    up <- 23.44 * (pi/180)
    around <- k * (pi/180)
    
     p <- 1
     x <- p * sin(up) * cos (around) 
     y <- p * sin(up) * sin (around)
     z <- p * cos(up)  
    
    xyz <- xyz.coords(x,y,z)
    polar_circle[1,k] <- as.numeric(xyz[1])
    polar_circle[2,k] <- as.numeric(xyz[2])
    polar_circle[3,k] <- as.numeric(xyz[3])
    polar_circle[4,k] <- 1
    
    }

#rgl.clear()
rotate_arctic <- 90*(pi/180)
rotate_antarctic <- 270*(pi/180)
arctic_x_L_rotation_matrix <- matrix(c(1,0,0,0,0,cos(rotate_arctic),-sin(rotate_arctic),0,0,sin(rotate_arctic),cos(rotate_arctic),0,0,0,0,1),4,4)
antarctic_x_L_rotation_matrix <- matrix(c(1,0,0,0,0,cos(rotate_antarctic),-sin(rotate_antarctic),0,0,sin(rotate_antarctic),cos(rotate_antarctic),0,0,0,0,1),4,4)
artic_circle <- arctic_x_L_rotation_matrix %*% polar_circle
antartic_circle <- antarctic_x_L_rotation_matrix %*% polar_circle

after_rotation
##                 x         y
## point_1 -4.832217 -1.284436
## point_2  4.827041  1.303754
sqrt((single_match[1,2]- single_match[1,1])^2+(single_match[2,2]- single_match[2,1])^2)
## [1] 7.071082
sqrt((after_rotation[1,2]- after_rotation[1,1])^2+(after_rotation[2,2]- after_rotation[2,1])^2)
## [1] 5.00003

Then translate that line onto the table

new_coord <- after_rotation + matrix(c(50,50,50,50),2,2)

Repeat for an entire box of matches

Throwing_matches <- function(){
  single_match <-cbind(c(-5,5),c(0,0))
  angle <- runif(1, min=1, max=360)*pi/180
  xy <- as.matrix(single_match)
  
  # Rotation
  cos.angle <- cos(angle)
  sin.angle <- sin(angle)
  after_rotation <- xy %*% t(matrix(c(cos.angle, sin.angle, -sin.angle, 
        cos.angle), 2, 2))
  
  # Translation
  translation.x <- runif(1, min=10, max=90)
  translation.y <- runif(1, min=10, max=90)
  translation <- matrix(c(translation.x,translation.x,translation.y,translation.y), 2, 2)
  new_coord <- after_rotation + translation
  
  
  return(new_coord)   
}

Throwing_matches()
##          [,1]     [,2]
## [1,] 23.69321 75.71892
## [2,] 25.25690 65.84194
for(i in 1:1000){
Throwing_matches()
}

surface_line <- function(line_height = 0){matrix(c(0,100,line_height,line_height), 2, 2, byrow=FALSE)}

  surface_line(20)
##      [,1] [,2]
## [1,]    0   20
## [2,]  100   20

#install.packages("retistruct")
library(retistruct)
one_match <- Throwing_matches()
one_surface_line <- surface_line(30)
line.line.intersection(one_match[1,], one_match[2,], one_surface_line[1,], one_surface_line[2,], interior.only = TRUE)
## [1] 14.58144 30.00000

Make a list of surface lines and a list of matches

list_of_matches <- array(NA,dim = c(2,2,number_of_matches_in_a_box))
list_of_surface_lines <- array(NA,dim = c(2,2,6))

for(i in 1:number_of_matches_in_a_box){
  list_of_matches[,,i] <- Throwing_matches()
}

for(i in 1:6){
  line_placement <- seq(0,100, by=20)
  list_of_surface_lines[,,i] <- surface_line(line_placement[i])
}

line_crossings <- matrix(NA,  dim(list_of_matches)[3] * dim(list_of_surface_lines)[3], 2)
counter <- 1
i <- 5
for(i in 1:dim(list_of_matches)[3]){
  for(j in 1:dim(list_of_surface_lines)[3]){
    
    line_crossings[counter,] <- line.line.intersection(list_of_matches[1,,i], list_of_matches[2,,i], list_of_surface_lines[1,,j], list_of_surface_lines[2,,j], interior.only = TRUE)
    
    counter <- counter + 1
  }
}

presence_of_line_crossing <- as.numeric(!is.na(line_crossings[,1]))

match_crosses_line <- length(which(!is.na(line_crossings[,1])==TRUE))
pi_estimate <- (2* match_stick_length * number_of_matches_in_a_box)/(distance_between_surface_lines * match_crosses_line)
pi_estimate
## [1] 3.194888